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ABSTRACT 

We use N-body simulations to investigate the radial dependence of the density, p, and velocity 
dispersion, er, in cold dark matter (CDM) halos. In particular, we explore how closely Q = 
p/o 3 , a surrogate measure of the phase-space density, follows a power law in radius. Our 
study extends earlier work by considering, in addition to spherically-averaged profiles, local 
Q-estimates for individual particles, Q t ; profiles based on the ellipsoidal radius dictated by 
the triaxial structure of the halo, Qi(r'); and by carefully removing substructures in order to 
focus on the profile of the smooth halo, Q s . The resulting Qf(r') profiles follow closely a 
power law near the center, but show a clear upturn from this trend near the virial radius, f2oo- 
The location and magnitude of the deviations are in excellent agreement with the predictions 
from Bertschinger's spherical secondary-infall similarity solution. In this model, Q oc r -1875 
in the inner, virialized regions, but departures from a power-law occur near r2oo because of 
the proximity of this radius to the location of the first shell crossing — the shock radius in the 
case of a collisional fluid. Particles there have not yet fully virialized, and so Q departs from 
the inner power-law profile. Our results imply that the power-law nature of Q profiles only 
applies to the inner regions and cannot be used to predict accurately the structure of CDM 
halos beyond their characteristic scale radius. 
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1 INTRODUCTION 



The study of the clustering of cold dark matter (CDM) on the scale 
of individual halos has progressed dramatically over the past cou- 
ple of decades due to the advent of powerful simulation techniques 
and ever faster computers. As a result, a number of basic prop- 
erties of the structure of CDM halos are generally agreed upon, 
even when many of these empirical findings lack a solid theoretical 
underpinning. One example is th e approximately "universal" mass 
profile of virialized CDM halos dNavarro et aljj 19961 . 1 19971 . here- 
after NFW). CDM halos are also strongly triaxial, with a preference 
for prolate shapes ( s ee, e.g., [Frenk et alii 19881 : Ijing & Sutcll2002l : 
lAllgood et al]|2006l : iHavashi et al.1 120071 , and refere nces therein) 
and have plentiful, albeit non-dominant, substructure ( Klyp in et al.l 
ll999l : lMooreetalJll999l) . 

The mass profile of CDM halos can be well approximated by 
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the simple law proposed by NFW, where the logarithmic slope of 
the spherically-averaged density profile follows the simple relation 
7 = rflogp/dlogr = —(1 + 3t)/(1 + x), with x — r/r_2 
being the radial coordinate expressed in units of a characteris- 
tic halo scale radius, r_2- More recent work shows evidence for 
small but systematic deviations from this simple law, and suggests 
that a third parameter may actually be requir ed to accurately de- 
scribe the mean mass profiles of CDM halos I Navarro et al.| 20041 : 



iMerritt et alj200ll2006l : [Gao et alj2008l : lHavashi & Whitei2008h . 
These authors argue that the density profile of ACDM halos steep- 
ens monotonically with radius, with no sign of converging to a 
central asymptotic power-law. The profiles are more accurately de- 
scribed by the "Einasto" formula for which 7 = — 2(r/r_2) a \ with 
a an adjustable "shape" parameter that can be tailored to provide 
an improved fit to a given halo density profile. 



As discussed by iNavarro et all J2008I) . this not only implies 
that CDM halos are not strictly self-similar, but also makes it diffi- 
cult to predict the asymptotic properties of CDM halo mass profiles. 
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For example, the Einasto and NFW formulae predict quite different 
asymptotic central behaviours: the Einasto profile has a true "core" 
with a well defined central density, whereas the central density of 
the NFW profile diverges l ike r . It seems clear from the lates t 
simulation results (see, e.g jNavarro et alj20 08; Stadel et alj2008h 
that the asymptotic slope is shallower than —1, but it is unclear 
whether the Einasto formula holds all the way to the centre and 
whether there is truly a well-defined central density for CDM ha- 
los aside from the ultimate upp er bound set by the Trem aine-Gunn 
phase-space density constraint dTremaine & GunrJI 19791) . The gen- 
tly curving nature of the Einasto profile makes it difficult to extrap- 
olate available simulations in order to predict the central properties 
of the halo with certainty. 

One alternative is suggested by the realization that, although 
the density, p(r), and velocity dispersion, a(r) are complex func- 
tions of radius, the quantity Q(r) = pi a 3 follows closely a sim- 
ple power-law, Q(r) oc r x , with roughly the sam e exponent, 

X ' 1.875, for all halos dTavlor & Navarro! 1200 ll) . Q(r) has 

the same dimensi ons as the phase-space density, /, but it is not a 
true m easure of it dAscasibar & B innev 2005; Sha rma & Steinmetzl 
l2006t) . We shall therefore refer to Q as a pseudo-phase-space den- 
sity, or as a surrogate measure of phase-space density. Despite 
this, Q relates two moments of / which occur often in equations 
that describe equilibrium systems, and therefore simple relations 
between them ar e extremely useful when constructing dynamical 
models (see, e.g. .[Austin et al ] 20051: IPehnen & McLaughlinll 20051 ; 
lBarnesetall2006l : lAscasibaretall2007h. 

The proposal of iTavlor & Navarro! d200ll hereafter, TN) 
has been confirmed in subsequent numerical work (see, e.g., 
Rasia e t al. 2004; Dehnen & McLaughlin 2 0051 ; iFaltenbacher et ail 
20071 : IVass et al. 2009; Wang & White 2009) and has been used in 



the literature to motivate dynamical models of dark halos. One in- 
triguing feature of the TN result is that the exponent of the power- 
law Q( r) profile is consisten t with that found in the similarity solu- 
tions of lBertschin ger ( 1985). Whether this is a mere coincidence or 
has a deeper meaning remains unclear. CDM halos form through a 
combination of smooth infall and the accretion of smaller progen- 
itors that are subsequently disrupted in the tidal field of the main 
halo. Bertschinger's similarity solution, on the other hand, follows 
the accretion of radial mass shells onto a point-mass perturber in 
an otherwise unperturbed Einstein-de Sitter universe. The solution 
assumes spheric al symmetry, allows only r adial motions, and is vi- 
olently unstable ( Vogelsbe rger et alj|2 009 ) . In spite of this, the ap- 
proximate power-law nature of Q(r) has been confirmed by the 
latest series of simul ations, which resol ve CDM halos with over 
one billion particles dNavarro et alfeooljh . 

The actual value of the exponent has also received attention. 
Although most simulations seem consistent with \ — —1-875, best 
fits often give slightly different values for y , typically in the range 
— 1.85 to —2 (but see I Schmidt et "al] |2008l for a differing view). 
Furthermore, there is an indication that y m ay depend on the main 
mode of mass accretion dWang & White 2009), and on the slope of 
the primordial power-spectrum dKnollmann et alj|2008l) . The cited 
work assumes that Q(r) is a power law, and then estimates x from 
simple fits to the spherically-averaged Q profiles. If Q profiles de- 
viate slightly but significantly from a power law, it could lead to 
a spread in the values of x< depending on, for example, the radial 
range of the fits or the characteristic mass of the halos considered. 

A further complication is introduced by the presence of sub- 
structure. Although not dominant in mass, because of their higher 
density and lower velocity dispersion, subhalos typi cally have 
much higher values of Q than the surrounding halo dArad et al.l 



l2004l : lDiemand et al.l2008l ; IVass et al.l2009l) . Together with the fact 
that CDM halos are in general triaxial, this hinders a proper defini- 
tion of the value of Q at given r, especially in the outer region s of 
a halo, where subhalos are most abundant dSpringel et al . 2008). 

We address these issues here using a series of high-resolution 
cosmological N-body simulations of the formation of individual 
CDM halos. In particular, we compute local estimates of Q at each 
particle position, Qi, and contrast the resulting profiles with those 
obtained using spherically-averaged estimates. The use of Qi al- 
lows us to carefully excise substructures and to focus our analysis 
on the pseudo-phase-space density profile of the smooth main halo. 

This paper is organized as follows. Sec.|2]provides a brief de- 
scription of our numerical simulations; Sec. [3] discusses our main 
findings and compares them with Bertschinger's similarity solu- 
tions. We conclude with a brief discussion of our main conclusions 
in Sec.g] 



2 THE SIMULATIONS 

We study the formation of 21 CDM halos selected from a 100 h~ 
Mpc-box cosmological simulation and resimulated at high res- 
olution in their full cosmological context. We provide below a 
brief summary of the numerical techniques, including the adopted 
cosmological parameters; the initial conditions setup; the simula- 
tion code; as well as the halo selection criteria and analysis tech- 
niques. More detailed information about our resimulation and anal- 
ysis techniques may be found in previous papers by our group 
(e.g..|PoweretalJl2003l : iNavarro et aT]|2004 ISpringel et al.ll2008l : 
iNavarro et al. 2008). 



2.1 Cosmological Parameters 

All our simulations adopt the currently favored ACDM cosmogony 
with the following parameters: Am = 0.25, Qa = 1 — Om = 0.75, 
erg = 0.9, n 3 — 1, and a Hubble constant Ho — H(z — 
0) = 100 ft km s" 1 Mpc -1 = 73 km s" 1 Mpc" 1 . These pa- 
rameter s are the same as tho se adopted for the Millennium Sim- 
ulation dSpringel et aT] |2005). and are consistent (within 2-sigma) 
with those derived from the WMAP 1- and 5-year data analyses 
(Sperg el et alj2003l : lKomatsu et alj2009l) an d with the recent clus- 
ter abundance analysis of Henry et all (2009). 



2.2 Halo Selection 

The 21 halos were selected from the same 900 3 -particle, 
100 h~ Mpc-box pa rent simulation used for the Aquarius Project 
dSpringel et alj|2008l) . These halos were subsequently resimulated 
at higher resolutio n using the technique described in detail by 
iPower et al.l d2003l) . We avoid halos that form in the periphery 
of much larger systems by imposing a mild isolation criterion so 
that no halo more massive than half the mass of the selected sys- 
tem lies within 1 h^ 1 Mpc at z = 0. This parent simulation 
was later also resimulated in its entirety at much higher resolu- 
tion; this is the Millennium -II Simulation recently analyzed by 
iBovlan-Kolchin et alj d200^) . 

Besides the six Aquarius halos, which were all selected to 
have virial[] masses of orders 10 12 h^Me , we have resimulated 



We define the virial mass of a halo, M200. as that contained within a 
sphere of mean density 200 times the critical density for closure, p cr it = 




Figure 1. Spherically averaged pseudo-phase-space density profiles for the 21 dark matter halos in our sample. The six level-2 Aquarius halos are shown at 
Z = (red dot-dashed lines), the other fifteen (solid blue) are shown at the most recent redshift when they pass the dynamical relaxation criteria (Sec. l2.4."3V 
Left panels correspond to Q(r) = p/er 3 ; panels on the right correspond to the "radial" Q r (r) = p/o" 3 . Radii are scaled to the scale radius, r_2, of each halo. 
Middle panels show residuals from the best r~ 1 - s75 power-law fit to the r collv < r < r_2 portion of the profiles. These best fits are also used to choose the 
vertical normalization of each profile in the upper panels, so as to minimize the halo-to-halo scatter in the inner profiles. Bottom panels are analogous to the 
middle ones, but for power-law fits over the whole range r conv < r < r200, with free-floating exponent, r x . Values of \ and \ r for each halo are listed in 
Table[2] 



a further set of 15 halos in order to span the mass range ~ 10 12 to a 
few times 10 14 ft -1 Mq. These 15 simulations have typically a few 
million particles within the virial radius at z — and are of lower 
numerical resolution than the level-2 Aquarius halos. Combining 
these two datasets allows us to assess the sensitivity of our results 
to numerical resolution. Table [T] lists the main properties of each 
halo in our sample. 



2.3 The Code 

All simulations were run with either the publicly available 
GADGET-2 code JSpringelll2005|) or its latest version, GADGET- 
3, which was developed for the Aquarius Project. Softening lengths 
are ch osen according to the "optimal" prescription of IPower et al.l 
d2003h . Pairwise interactions are fully Newtonian for separations 
exceeding the spline-lengthscale h s . Table Q] quotes the equiva- 
lent Plummer softening, ec = h 3 /2.8, for each resimulated halo. 
Throughout the simulations, the softening length is kept fixed in 
comoving coordinates. 



2.4 Analysis 

2.4.1 Spherically-averaged Q profiles 

In order to compute the spherically-averaged pseudo-phase-space 
density profile of each halo we first identify the halo center with 

3H 2 /8nG. The virial mass defines implicitly the virial radius, r2oo. and 
virial velocity, V200 = (GM200 A200) 1 , of a halo, respectively. 



the location of the particle having the minimum potential energy. 
Then we compute Q(r) in iVbm spherical shells equally spaced in 
log 10 r in the range rco nv ^ r ^ r 200- H ere r conv is the conver- 
gence radius defined bv lPoweret all 42003^ where circular veloci- 
ties converge to better than 10% (see also lNavarro et al.l2008l) . For 
each spherical shell (radial bin), we estimate Q(r) = p/cr 3 , where 
p is just the mass of the shell divided by its volume and a 2 is twice 
the specific kinetic energy in the shell. We also compute a "radial" 
Q estimate, Q r (r) = p/cr 3 in an analogous way, although instead 
of the total kinetic energy we use only the kinetic energy in radial 
motions to estimate ay. 



2.4.2 Local Q profiles 

A different estimate of the pseudo-phase-space density may be ob- 
tained, for each shell, by considering "local" estimates of Q at the 
position of each particle. We shall call this Qi — pi/crf, and use 
JVngb nearest neighbours in order to compute the local density and 
velocity dispersion. Density estimates at the location of the i th par- 
ticle are computed as 



A, hi), 



(1) 



where n,j = Vi — rj, and W(r, h) is the smoothing kernel often 
adopted in Smoothed Particle Hydrodynamics (SPH) simulations: 



W(r, h) 




< - < - 

U h ^ 2' 
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Figure 2. Radial profiles of the velocity anisotropy parameter, (8 = 1 — 
a 2 /2a 2 ,, for all halos in our sample. Line types are as in Fig.[T] Radii are 
normalized to the scale radius of each halo. Note the non-monotonic radial 
dependence of (3: halos are nearly isotropic near the center, radially biased 
near r_2, but approximately isotro pic again in the outskirts. The anisotropy 
expected from the ,8(7) relation of Hansen & Moore 12006) for an Einasto 
profile with a = 0.174 is also shown. 

The smoothing length, hi , of each particle is defined implicitly by 
the smallest volume that contains N ng b nearest neighbours: 

h ^{ N ^y 3 - (3) 

We use, as default, iVngb = 48 for our lower resolution runs and 
64 for the level-2 Aquarius halos. 

Given iV n gb, the local velocity dispersion for particle i is given 
by of = v 2 — v 2 , where the unweighted averages are computed 
over all iV ng b neighbours. 

2.4.3 Relaxation criteria 

In order to minimize the effect of transient, rapidly-evolving evo- 
lutionary stages, such as ongoing mergers, we impose (when ex- 
plicitly sta t ed) re laxation criteria similar to those introduced by 
iNeto et al.l fe007t) . These include restrictions on the fraction of 
the virial mass in self-bound substructures, / su b = M su b(r < 
7"2oo)/A^20o < 0.07; on the offset between the center of mass of 
the halo and its true centre (as defined by the particle with mini- 
mum potential energy), d g — \vcm — r ccn \/r2oo < 0.05, and on 
the virial ratio of kinetic to potential energies, 2i\7|<I>| < 1.3. 

In practice, when a halo does not satisfy the relaxation crite- 
ria at z — we track its main progenitor back in time until we 
find the first snapshot when it does. This typically occurs at red- 
shifts less than ~ 0.2, but in one case we had to go back in time 
until z ~ 0.8 in order to find a suitably "relaxed" configuration. 
In what follows, we shall consider relaxed configurations only for 
the lower-resolution halos but take the z — configuration for the 
Aquarius halos. As we show below, the results are similar in the 



Figure 3. "Local" estimates of the pseudo-phase-space density, Qi, as a 
function of distance from the halo centre for all halo particles in one of our 
simulations. Different colors correspond to various particle subsamples. (i) 
Blue denotes particles in self-bound substructures, as identified by SUB- 
FIND, (ii) Cyan indicates particles not bound to any substructure but with 
higher-than-average Qi for their location in the halo. These are particles in 
recently- stripped tidal streams which, although assigned to the main halo by 
SUBFIND, have yet to phase mix within the main halo i Macieie wski et all 
2009). (iii) Red dots indicate particles with Qi < Q C ut(r), which we de- 
fine as belonging to the relaxed main halo (see Fig. [4}- Also plotted are the 
spherically-averaged Q(r) profile (solid blue line), and the best-fit r -1 - 875 
power law (dot-dashed line). Vertical and horizontal bands show, respec- 
tively, the particles selected for Figs.|4]and \5\ 

two cases, which means that our conclusions are not particularly 
sensitive to our requirement of dynamical equilibrium. 

The properties of each halo in our sample at z — are listed 
in Table [T] Here we list the virial mass, M200, the virial radius, 
^200, the number of particles (./V200) within V2oo, as well as the 
gravitational softening, ec, and the convergence radius, r conv . The 
peak of the circular velocity curve is also specified by r max and 

Vmax • 



3 PSEUDO-PHASE-SPACE DENSITY PROFILES 
3.1 Spherically-averaged Q profiles 

Figure [T] shows the spherically-averaged Q(r) profiles for all ha- 
los in our sample, together with residuals from various best-fits. 
Panels on the left and on the right correspond to Q(r) and Q r (r), 
respectively. The plotted profiles extend from the convergence ra- 
dius, r conv , to the virial radius, r2oo- Middle panels show residuals 
from best fits to the region inside r_2 with an r -1 ' 875 power law. 
All profiles are normalized to the scale radius, r_2, and vertically 
according to the power-law best fit. Bottom panels show residuals 
from fits to the r conv < r < r2oo profile with a power law with 
free-floating exponent, Q oc r x . 

A few things are worth noting in this figure. The first is how 
closely both the Q(r) and Q r (r) profiles follow simple power 
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Figure 4. Distribution of local phase-space densities, Qi, for particles in 
the thin (shaded) spherical shell near the virial radius of the halo shown in 
Fig. [3] Substructures identified by SUBFIND are shown in blue, the main 
SUBFIND halo in cyan. Various thin lines illustrate the effect of varying the 
number of neighbours in the pi and Oi estimates, as labelled in the legend. 
As discussed in the text, a simple cut Qi < Q C ut identifies unequivocally 
all well-mixed particles in the main halo; Q C ut is shown by the vertical 
dotted line. Down-pointing arrows indicate the median of the distribution 
of each set of particles. Because there are few particles in the high-Q ; tail, 
the median Qi of the main SUBFIND halo and that of particles with Qi < 
Q cu t are nearly identical. We adopt the latter as the characteristic pseudo- 
phase-space density of the smooth halo at each radius. 



laws, from the innermost resolved radius out to r2oo- In the case 
of Q(r), even when the exponent of the fit is fixed at \ — —1-875, 
which means that a single free parameter (the vertical normaliza- 
tion) is allowed, residuals from best fits do not exceed ~ 30% any- 
where within the virial radius. This power-law behaviour holds for 
roughly three decades in r and six decades in Q. 

Although the best-fit x differs from —1.875 (see Table|2]for 
actual values), the residuals decrease only very slightly when al- 
lowing x t° float freely. Defining a figure-of-merit function as 

^ 2 = tt- ^ (InQ-lnQflt) 2 , (4) 

we find that %[>, averaged over all halos and fit over the range 
fconv < r < V2oo, varies only from (ip) — 0.102 when fixing 
X = —1.875 to (tp) — 0.085 when allowing x to be a free param- 
eter. 

The "tilt" in the Q r -residuals shown in the middle-right panel 
of Figure fJJ indicates that the Xr exponent that fits best the Q r (r) 
profiles is slightly more negative than —1.875. Overall, however, 
Q r (r) may also be approximated with a power law with Xr ~ 
1.92, as may be seen in the bottom right panel of the same figure. 
The average xp for power-law fits with variable Xr is 0.142, which 
means that Q r (r) deviates more than Q(r) from a simple power 
law, but only slightly so. 

The reason why Q r (r) deviates more from a simple power 



law than Q(r) may be traced to the complex behaviour of the 
anisotropy parameter, /3(r) = 1 — erf /2of . All halos in our sam- 
ple are almost isotropic near the center, radially anisotropic further 
out, but nearly isotropic again close to the virial radius. Because 
of this complex radial behaviour, Q(r) and Q r (j) cannot both be 
simultaneously well fitted by a power law with the same exponent. 

Although power laws provide excellent fits to the pseudo- 
phase-space density profiles, further scrutiny of the middle and 
lower panels of Fig. fJJ reveals a well defined trend in the residuals 
of most halos, which tend to "curve up" slightly but significantly 
in the outer regions and, to a lesser extent, in the innermost regions 
as well. The latter deviations are best appreciated in the Aquarius 
halos, which have much better resolution than the rest. 

These results seem to apply both to the Aquarius halos and to 
the dynamically-relaxed halo sample, which suggests that our con- 
clusions are not crucially dependent on the adoption of the particu- 
lar relaxation criteria we used to select the sample. The above-noted 
trend in the residuals means that the exponent, \, derived from r ~ x 
power-law fits will depend on the radial range adopted for the fit. 

Given the desirable properties of a simple power law, it is 
worth investigating whether the deviations from a simple r x be- 
haviour might be due to the presence of substructure or to the 
aspherical nature of halo structure. We explore these possibilities 
next. 

3.2 Local Q profiles 

Fig. f5] shows Qi, the local estimate of Q at the location of each 
particle in the halo, as a function of the distance from the halo cen- 
tre. Because substructures are overdense and have lower velocity 
dispersion than their immediate surroundings, they show up promi- 
nently in this plot as particles with very high Qi at a given radius. 
This is confirmed by the color coding adopted in the figure: par- 
ticles in d ark blue are those a ssociated by the substructure-finder 
SUBFIND dSpringel et alj200ll) to self-bound subhalos that survive 
within the main halo. Clearly, because their pseudo-phase-space 
density is so distinct from the main halo, substructures have the 
potential to bias estimates of Q(r), especially in the outer regions, 
where subhalos are more prevalent. 

Although it would be simple enough to remove the self-bound 
structures from Q(r) profiles, the cyan dots in Fig. |3] illustrate a 
second, related problem. These are particles that SUBFIND asso- 
ciates with the main halo, but which clearly have deviant Qi values 
relative to the surroundin g average. As discussed, for example, by 
iMacieiewski et al] ( l2009l) . these are particles recently stripped from 
substructures; although now unbound to any subhalo, they have yet 
to phase mix fully with the underlying main halo. 

Fig. E] shows the Qi distribution of all particles in the thin 
spherical shell near the virial radius of the halo shown by the verti- 
cal shaded region in Fig.[3] This shows the wide range in Qi (almost 
8 decades) spanned by particles at a fixed distance from the halo 
centre. Despite this, the figure also shows that the characteristic Q 
at that distance is well defined, since most particles have Qi values 
close to the peak on the left of the distribution. (Note the logarith- 
mic scale in the j/-axis.) The tail to the right of this peak contains 
self-bound substructures (in blue, as identified by SUBFIND) and re- 
cently stripped material, which, as mentioned above, are assigned 
to the main halo by SUBFIND despite their higher-than-average Qi 
values. 

The shape of the distribution suggests that imposing a simple 
criterion, e.g., Qi < Qcut, should be sufficient to identify unequiv- 
ocally well-mixed particles belonging to the main smooth halo; a 
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Figure 5. Orthogonal projections of particles in a shell of roughly constant pseudo-phase-space density (i.e., those in the horizontal shaded band in Figure|5] 
excluding substructures). For clarity, we plot in each panel only particles in a thin slice perpendicular to the viewing axis. Red points show the original 
positions; these delineate a nearly prolate ellipsoid, which has been rotated so that its principal axes coincide with the coordinate axes of the plot. Green points 
show the projected particle positions after correcting for triaxiality by the method outlined in the text. 



plausible choice for Q cut is shown by the vertical dotted line in 
Fig. [4] Down-pointing arrows indicate the median of the distribu- 
tion of each set of particles. Note that, because a small fraction of 
the particles are in the high-Qi tail, the median Qi of the main 
SUBFIND halo and that of particles with Qi < Q cu t are nearly 
identical. 

We therefore adopt the simple Qi < Q C ut prescription to de- 
fine the main "smooth" halo. Once Q C ut is chosen at some radius, 
we may use the approximate power-law behaviour of Q to scale it 
to any other radius by Q cu t(r) oc r -1,875 . Particles shown in red 
in Fig.[3]are those assigned to the smooth halo by this criterion. At 
each radius, we shall adopt the median Qi of these particles in or- 
der to construct the pseudo-phase-space density profile of the main 
smooth halo. We shall refer to this profile as Ql(r). 

We note that this definition is insensitive to the number of 
nearest neighbors (-/V ng b) adopted to compute Qc the various lines 
in Fig.|4]illustrate the results for iV n gb = 48 particles (our default 
value), as well as for 1000 (long-dashed), 500 (dotted), and 100 
(dot-dashed), respectively. 



3.3 Correction for triaxiality 

Dark halos are not spherically symmetric. Because of this, a shell of 
particles at constant distance from the halo center will have a wide 
Qi distribution, even if one subtracts substructure as specified in the 
previous subsection. Iso-Q surfaces track fairly well the isodensity 
contours of the main halo, and follow closely three-dimensional el- 
lipsoidal surfaces. Fig.[5]shows three orthogonal projections of par- 
ticles with similar values of Qi, selected from those falling in the 
horizontal band highlighted in Fig. [3] Only particles in the smooth 
main halo are plotted here. The original particle positions (in red) 
in a thin slice perpendicular to the line of sight are seen to trace 
a nearly prolate ellipsoid, which for convenience has been rotated 
so that its principal axes coincide with the coordinate axes of the 
projection. 

Given that the "iso-Q" surfaces are well approximated by el- 
lipsoids, we may use the eigenvalues of the diagonalized inertia 
tensor to compute an elliptical radius, r', for each Qi, to define 
an "ellipsoidal" Qi(r') profile that may be contrasted directly with 
Qi (r). In practice, we slice the smooth main halo in narrow bins in 
Qi; compute the axis lengths a, b, and c; and use them to reassign 
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Figure 6. Local pseudo-phase-space density as a function of radius for the 
same halo shown in Fig. [3] Colored dots are as in Fig. |3] Red and green 
dots indicate smooth main halo particles before and after correcting for halo 
shape. The solid curve is the spherically-averaged profile, and is compared 
with the curves tracing the median Qi for (i) all particles (blue dot-dashed), 
(ii) all particles in the smooth main halo (red dashed), and (iii) smooth main 
halo particles corrected for triaxiality (dashed). The dotted line shows a 
power-law, r~ 1 - 875 . 



an ellipsoidal radius r to each particle in the smooth main halo. 
We compute r' as 



(5) 



after normalizing a, b, and c so that abc = 1 for all shells. This 
choice preserves the typical distance to the halo centre of particles 
in a given Qi shell. The result of the procedure is illustrated in 
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Figur e 7. Spherically-ave raged (left) and median local (right) pseudo-phase-space density profiles for the convergence series of the Aq-A halo jSpringel et alj 
2008; Navarro et al. 2008). Different colors correspond to runs with different resolution, as indicated in the legend. Local profiles refer to the median Qi of the 
smooth main halo as a function of distance to the halo centre. The dotted line is the power-law Q <x r -1 875 , scaled to match Aq-A-2 at r ^ r_2 = 11.15 /i — 1 
kpc. Bottom panels show the residual with respect to this power-law. 



Fig. [5] where the green dots are the same particles as those shown 
in red, but the coordinate axes are now x' = x/a, y' — y/b and 
z' = z/c. 

The green dots in Fig. [6] delineate the Ql(r') profile for the 
smooth main halo. The lines in Fig. [6] show the variations in the 
pseudo-phase-space density profile induced by the various alterna- 
tive ways of estimating distances and Q that we have discussed so 
far. 

Although the profiles change appreciably, they all show the 
same upturn in the outer regions, relative to a simple j-~ 1,875 
power-law fit, noted when discussing the spherically-averaged 
Q(r) profiles (Sec. 13. lb . The upturn is indeed more pronounced 
when using local-Q estimates, because (i) local densities are sen- 
sitive to inhomogeneities and are therefore higher than the spher- 
ical average, and (ii) because local a estimates are lower than the 
spherical average, since they do not include the bulk motion of sub- 
halos and recently stripped material. Interestingly, Fig. [6] shows 
that, after correcting for triaxiality, the spherically-averaged profile 
is indistinguishable from the Qi(r') profile. We conclude that the 
upturn is not caused by the presence of substructures in the outer 
regions nor by departures from spherical symmetry. We discuss the 
interpretation of this robust feature of the Q profiles next. 



3.4 Numerical convergence 

Before considering the meaning of the departures of Q profiles 
from simple power laws we should check explicitly that our con- 
clusions are not affected by the numerical resolution of the simu- 
lations. We show this in Fig. [7] where we compare the Q(r) and 
Ql (r) profiles of one of the Aquarius halos at four different reso- 
lutions. The highest, Aq-A-2, has more than 100 million particles 
within the virial radius; the lowest, Aq-A-5, about 600 thousand. 
The non-Aquarius halos in the series analyzed in this paper have 
numerical resolution comparable to Aq-A-4. Fig. UJ shows con- 



vincingly that our results are unlikely to be adversely affected by 
numerical resolution. Both the spherically-averaged and the local 
pseudo-phase-space density profiles are very well reproduced at all 
radii, down to the inermost resolved radius, r conv , of each run. 



3.5 Comparison with Bertschinger's similarity solution 

The local pseudo-phase-space density profiles for the smooth main 
halo of all our systems is shown in Fig. [8] The profiles are shown as 
a function of the elliptical radius (eq.[5j, scaled to the scale radius of 
each halo, r_2- All profiles have been normalized vertically so that 
they coincide at r' = r_2. The dotted line shows a Q oc j- -1 - 875 
power law, also normalized at r_2- The bottom panel shows resid- 
uals relative to the r - 1,875 power law. 

Fig. [8] shows the main result of our analysis. The pseudo- 
phase-space density profiles of our simulated halos clearly devi- 
ate from a simple power law in the outer regions. This devia- 
tion is actually pred i cted b y the secondary infall similarity solu- 
tion of iBertschingej {T985). Indeed, in the similarity solution the 
Q oc r - 1,875 behaviour occurs only asymptotically in the inner 
regions of the halo. In the outer regions, and more precisely, near 
the location of the "shock radius", r s h oc k (for collisional fluids), or, 
equivalently, of the first infall caustic (for collisionless fluids), the 
pseudo-phase-space density profile (or entropy profile in the case of 
a collisional fluid) shows a clear upturn from the inner asymptotic 
power-law behaviour. 

The reason for this upturn is that the fluid is not fully viri- 
alized near r a hock, since this radius marks the transition between 
mass shells that are infalling for the first time and those that have 
already crossed (or shocked) material that collapsed earlier. For ex- 
ample, in the case of collisionless fluids, a mass shell must cross 
the center and complete roughly 2-3 full oscillations before settling 
onto a periodic orbit of constant apocenter. In the case of a colli- 
sional fluid, a newly shocked shell drifts inward from the radius at 
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Figure 8. Local estimates of the pseudo-phase-space density as a function 
of the ellipsoidal radius r' (eq.|5J. Each curve traces the median Qf(r') 
for smooth main halo particles, computed in r' bins of equal logarithmic 
width. Radii are scaled to the value of r_2 of each halo. The dotted line 
shows the Q oc r — 1-875 p 0wer l aWj whereas the thick dashed black curve 
shows the actual similarity solution computed by Bertschinger 1 1985). The 
outermost point of this solution, indicated by the downward pointing arrow, 
corresponds to r s hock> th e "shock radius" (for collisional fluid infall) or 
to the position of the first caustic (for collisionless fluids). Radii for the 
similarity solution assume rg^ock = 8r_2 All profiles are normalized 
at r' = r—2- Note that halo profiles deviate from the asymptotic inner 
r — 1.875 p OW£r j aw m me same wa y as the similarity solution. 



which it was shocked before reaching hydrostatic equilibrium (see 
Fig. 4 of iBertsch inger 1985). As a result, Q(r) (or the "entropy" 
in the case of a collisional fluid) shows a characteristic upturn at 

Tahock like the one sh own in Fig. [8] 

As discussed by IWhite et al] ( fl993h , the first caustic/shock, 
which occurs at about a third of the current turnaround radius, lies 
close to the virial radius, as defined here. Giv en that our halos have 
concentrations, c = r2oo /v~2, of order 7-10 dNeto et aT1l2007t) . we 
would then expect the upturn in the Q profiles to occur roughly at 
~ 7-10 r_ 2 . 

The dashed black curves in Fig.[8]show the similarity solution, 
assuming r s h oc k = 8 r_2, and normalized vertically at r = r_2 to 
coincide with the simulated halo profiles. It is clear from Fig.[8]that 
the similarity solution is in excellent agreement with the Q profiles 
of our simulated halos. This suggests that the upturn in the outer Q 
profiles just reflects the fact that the regions around the virial radius 
have not yet fully virialized. 



4 SUMMARY 

We have used a set of 21 high-resolution cosmological N-body 
simulations to investigate the pseudo-phase-space density profiles, 
Q(r) = p/cr 3 , of cold dark matter halos. In particular, we concen- 
trate our analysis on the radial dependence of Q for particles in the 



smooth main halo component, after carefully removing substruc- 
tures and correcting for halo shape. 

Our main result is that although Q(r) is remarkably well- 
approximated by a power law, a slight but systematic upturn from 
the power-law profile is clearly seen in the outer regions of all our 
simulated halos. Both the exponent of the power law, as well as the 
upturn in the outer regions are consistent w ith the particu l ar sec - 
ondary infall similarity solutions derived by Bertschi nger] d 1985b - 
In these solutions, the power-law inner region corresponds to the 
virialized region of the halo, whereas the upturn in the outer re- 
gions coincides with the location of the shock/first caustic of the 
system or, roughly speaking, with the virial boundary of a halo. Al- 
though Bertschinger's is just one particular secondary-infall simi- 
larity solution, valid in the simple case of a point-mass perturber 
in an Einstein-de Sitter universe, the upturn in Q marking the virial 
boundary of the system is expected to be a general feature of such 
solutions. 

Although the outer upturn may be robust, the fact that 
the exponent of the inner power law is consistent with 
Be rtschinger's (y ~ —1 . 875) i s somewhat surprising. As shown 
bv lFillmore & Goldreichl (l984), the non-linear structure of halos 
formed through self-similar secondary infall depend on the scaling 
index that characterizes the mass dependence of the initial perturba- 
tion, SM /M oc A/~ E . In the case of Bertschinger's solution e = 1, 
but this is a poor approximation to the typical overdensity that seeds 
the collapse of galaxy-sized LCDM halos. The agreement between 
Bertschinger's Q(r) and the simulated profiles is therefore a non- 
trivial result whose significance remains unclear. 

The presence of the upturn in the Q(r) profiles casts doubts 
on work that attempts to construct dynamical equilibrium mod- 
els of CDM halos by assuming that the p ower-law behaviour 
of Q profiles applies t o all radii (see, e.g., lAustin et al.l I2005L 
Dehnen & McLaughlin' 120051) . The r x behaviour seems to hold 
only in the inner regions, and our results caution against fitting 
power laws to Q(r) over a radial range that extends outside the 
scale radius. 

Because the radius where the upturn becomes noticeable 
marks the transition to the region where virial equilibrium no longer 
holds, this radius, expressed as a fraction of the virial radius, may 
vary systematically with halo mass, collapse time, or cosmological 
parameters. Indeed, the virial radius definition we adopt here does 
not depend on the formation history of each object, but the bound- 
ary of the region where virial equilibrium might hold does. For ex- 
ample, the shock radius of an early collapsing halo that has accreted 
little mass in the recent past might occur further away (in units of 
i"20o) than another system of similar mass that has assembled a con- 
siderable fraction of its present-day mass more recently. It may be 
worth checking whether the trends of \ with po wer spectrum and 
mode of assembly noted in t he literature (see, e.g., Kn ollmann et al.l 
2008: Wang & Whitel2 009) might be explained by this process. 

Finally, we note that our highest resolution halos (those from 
the Aquarius project) also present evidence of departures from a 
simple power-law Q profile near the innermost resolved radius. At 
the moment it is unclear whether this indicates that the power-law 
behaviour of Q(r) is just an approximation that breaks down inside 
some characteristic radius, or whether estimates of Q(r) at those 
radii might be affected by numerical uncertainties. For example, it 
is easy to demonstrate that an isotropic halo whose mass profile 
follows strictly an Einasto law cannot have a pow er-law Q(x) law 
that extends all the way to the centre (see, e.g., iMa et al. 2009). 
The radial dependence of these Q(r) profiles may be more accu- 
rately represented by an Einasto-like form where the power-law in- 
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dex changes gradually but smoothly with radius. We plan to address 
this and other pertinent issues in forthcoming work. 
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Table 1. We list, for each halo in our sample, a label, the Plummer-equivalent gravitational softening and convergence radius, as well as the virial radius at 
z = 0. The structural parameters V ma x and r max identify the peak of the circular velocity curve for each halo. Finally, M200 is the virial mass and -/V200 
is the number of particles within the virial radius. The z Ici column specifies, for non-Aquarius halos, the redshift of the most recent snapshot when the halo 
satisfies the relaxation criteria introduced in Sec. 12.4.31 



Halo 








f conv 


T200 


T max 


Vmax 


M200 


N200 








[kpc/h] 




[kpc/h] 


[kpc/h] 


[kpc/h] 


[km/s] 


[lO^Me/h] 


[106] 




Aq-A-2 


4 


8x10" 


2 


0.250 


179.5 


20.5 


208.5 


134.5 


134.5 




Aq-B-2 


4 


8x10" 


2 


0.219 


137.0 


29.3 


157.7 


59.8 


127.1 




Aq-C-2 


4 


8x10" 


2 


0.248 


177.3 


23.7 


222.4 


129.5 


126.8 




Aq-D-2 


4 


8x 10~ 


2 


0.281 


177.3 


39.5 


203.2 


129.5 


127.0 




Aq-E-2 


4 


8x 10~ 


2 


0.223 


155.0 


40.5 


179.0 


86.5 


123.6 




Aq-F-2 


4 


8x10" 


2 


0.209 


152.7 


31.2 


169.1 


82.8 


167.5 




hi 




0.39 




1.090 


134.4 


44.0 


151.9 


56.4 


2.04 


0.198 


h2 




0.25 




0.852 


144.6 


35.1 


159.9 


70.3 


4.91 


0.000 


h3 




0.38 




1.063 


154.1 


33.7 


178.6 


85.0 


2.60 


0.049 


h4 




0.31 




0.899 


154.7 


30.3 


175.8 


86.1 


4.17 


0.876 


h5 




0.24 




0.797 


156.1 


32.0 


174.8 


88.5 


6.10 


0.000 


h6 




0.26 




0.829 


158.0 


47.3 


171.6 


91.7 


6.00 


0.000 


hi 




0.35 




1.001 


158.3 


37.4 


184.8 


92.2 


3.25 


0.062 


h8 




0.39 




1.148 


175.6 


39.1 


203.5 


125.9 


3.30 


0.350 


h9 




0.45 




1.351 


177.8 


45.1 


200.2 


130.8 


2.48 


0.000 


hlO 




0.33 




0.944 


183.7 


21.3 


209.2 


144.1 


4.79 


0.350 


hi I 




1.06 




3.239 


275.5 


188.2 


285.6 


486.4 


1.08 


0.309 


hl2 




1.39 




4.006 


391.0 


114.2 


425.6 


1389.5 


1.26 


0.140 


hl3 




1.61 




4.418 


396.1 


86.9 


422.9 


1445.2 


0.98 


0.140 


hl4 




2.08 




6.886 


856.0 


263.6 


888.2 


14581.8 


2.62 


0.030 


hl5 




3.65 




11.595 


981.6 


444.8 


1043.0 


21993.1 


1.16 


0.094 



Table 2. Values of the exponent \ obtained from r x fits to various Q profiles, as indicated in the legend. Xr refers to fits to the "radial" Q r profiles. The 
average \ f° r a 'l halos is listed in the next-to-last row, together with its standard deviation . The average figure-of-merit for all halos, {i/>min)> lS listed in the 
last row of the table. 
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-1.846 


-1.792 


hlO 


-1.975 


-2.042 


-1.871 


-1.991 


-1.972 


-1.754 


hll 


-1.896 


-2.008 


-1.843 


-1.845 


-1.903 


-1.756 


hl2 


-1.879 


-1.949 


-1.875 


-2.023 


-1.864 


-1.758 


hl3 


-1.843 


-1.955 


-1.847 


-1.882 


-1.768 


-1.724 


hl4 


-1.810 


-1.902 


-1.817 


-1.886 


-1.801 


-1.732 


hl5 


-1.847 


-1.966 


-1.835 


-1.919 


-1.765 


-1.648 


(x) 


-1.893 ±0.048 


-1.963 ±0.050 


-1.858 ±0.031 


-1.919 ±0.047 


-1.876 ±0.061 


-1.777 ±0.052 


<^min> 


4.3(±1.3) x 10~ 2 


5.38(±0.16) x 10- 1 


8.46(±2.1) x 10~ 2 


1.42(±0.50) x 10- 1 


4.23(±1.79) x 10~ 2 


1.18(±0.43) x 10- 1 



